##### Compute GMM Objective Value
	
	function compute_gmm_objective(theta)
	
		theta_CLMS = theta[1:length(claims_moments_covariates)];
		avg_claims_moment_values = transpose(X_claims_data .* (ac_weights ./ sum(ac_weights))) * 
			(exp.(y_claims_data) .- exp.(X_claims_data * theta_CLMS));
		foc_moment_values = compute_foc_moments(theta,0); 
		moment_values = cat(dims=1,avg_claims_moment_values,foc_moment_values);
	
		return(transpose(moment_values) * W * moment_values)
	end
	
	function compute_optimal_weight_matrix(theta)
	
		theta_CLMS = theta[1:length(claims_moments_covariates)];
		
		W = zeros(num_moment_conditions,num_moment_conditions);
		foc_moment_contribution_values = 
			#foc_moment_contributions!(cat(dims=1,theta_CLMS,theta_RS_initial),zeros(length(plan_pmt_names)));
			foc_moment_contributions!(theta_CLMS,zeros(length(plan_pmt_names)));
		plan_pmt_demand = pmt[:,findall(pmt_fields .== "real_demand")[1]];
		
		for plan_to_check in plans_pmt[plans_pmt[!,:year] .<= year,:pmt_name]
			j = findall(plans_pmt[!,:pmt_name] .== plan_to_check)[1];
			moment_contribution = zeros(num_moment_conditions);
			moment_contribution[1:length(theta_CLMS)] = get_claims_moment_contribution_values(j,theta_CLMS);
			moment_contribution[findall(foc_moment_cond_names .== plan_foc_moment_cond_names[j])[1] + length(theta_CLMS)] = 
				foc_moment_contribution_values[j];
			W[1:num_moment_conditions,1:num_moment_conditions] += 
				(moment_cont_weights[j] ./ sum(moment_cont_weights)) .*
				#(plan_pmt_demand[j] ./ 
				#	sum(plan_pmt_demand[findall(plan_foc_moment_cond_names .== plan_foc_moment_cond_names[j])])) .* 
				(moment_contribution * transpose(moment_contribution));
		end
		
		return(inv(W))
	end	
	
	function compute_standard_errors(theta)
	
		theta_CLMS = theta[1:length(claims_moments_covariates)];
		
		# Nice discussion in Greene on pp.538-539
	
		# Compute Weight_matrix (M * M, where M is number of moment conditions)
		Weight_matrix = compute_optimal_weight_matrix(theta);
		
		# Compute G (M * K, where K is the number of parameters)
			# Need derivative of each moment with respect to each parameter
		
		G = zeros(num_moment_conditions,length(theta));
		
		# Add AC moments
		for m in 1:length(theta_CLMS)
			compute_ac_moment_value(theta_CLMS) = compute_avg_claims_moment_value(theta_CLMS,m);
			G[m,1:length(theta_CLMS)] = Calculus.gradient(compute_ac_moment_value,theta_CLMS);
		end
		
		# Add FOC moments
		for m in 1:length(foc_moment_cond_names)
			compute_foc_moment_value(theta) = compute_foc_moments(theta,m); 
			G[(m+length(theta_CLMS)),:] = Calculus.gradient(compute_foc_moment_value,theta);
		end
		
		vcov_matrix = inv(transpose(G) * Weight_matrix * G);
		#standard_errors = sqrt.(1/length(plan_names) *  diag(vcov_matrix));
		standard_errors = sqrt.(diag(vcov_matrix));
		
		return(standard_errors)
	end
	
##### Moment Conditions

	function foc_moment_contributions!(theta,foc_moment_contribution_values)
			
		# Calculate partial_firmclaims_partial_price based on candidate theta
		average_claims = exp.(X_CLMS * theta);
		ac_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * average_claims;
		demand_partial_ac_partial_price = theta[findall(claims_moments_covariates .== "log_risk_score")[1]] .*
			(transpose(rs_partial_matrix .* own_plan_year_matrix) * (average_claims .*
				plan_demands ./ plan_risk_scores));
		partial_firmclaims_partial_price_theta = (1 ./ partial_plan_demand_partial_price) .*
			ac_partial_ownplans_partial_price .+ demand_partial_ac_partial_price;
		
		foc_moment_contribution_values[1:length(plan_pmt_names)] = 
			partial_firmclaims_partial_price_focs .- partial_firmclaims_partial_price_theta;
	end

	function get_claims_moment_contribution_values(j,theta_CLMS)
	
		plan_name = plans_pmt[j,:plan_name];
		plan_year = plans_pmt[j,:year];
		pmt_name = plans_pmt[j,:pmt_name];
		pt_name_zero = string(plan_name,plan_year,0);
		in_avg_claims = findall(in(moments_data[!,:pmt_name]).(plans_pmt[!,:pmt_name]));
		plan_pmt_demand = pmt[:,findall(pmt_fields .== "real_demand")[1]];
		
		if pmt_name in moments_data[!,:pmt_name]
			data_index = findall(moments_data[moments_data[!,:year] .<= year,:pmt_name] .== pmt_name)[1];
			moment_contribution_values = 
				exp.(X_claims_data[data_index,:]) .* 
					(exp.(y_claims_data[data_index]) .- exp.(transpose(X_claims_data[data_index,:]) * theta_CLMS));
		elseif pt_name_zero in moments_data[!,:pmt_name]
			data_index = findall(moments_data[moments_data[!,:year] .<= year,:pmt_name] .== pt_name_zero)[1];
			common_indices = setdiff(intersect(findall(plans_pmt[!,:plan_name] .== plan_name),
					findall(plans_pmt[!,:year] .== plan_year)),in_avg_claims);
			moment_contribution_values = (plan_pmt_demand[j]/sum(plan_pmt_demand[common_indices])) .*
				exp.(X_claims_data[data_index,:]) .* 
					(exp.(y_claims_data[data_index]) .- exp.(transpose(X_claims_data[data_index,:]) * theta_CLMS));
		else
			moment_contribution_values = zeros(length(theta_CLMS));
		end
		
		return(moment_contribution_values)
	end

	function compute_foc_moments(theta,m)
		
		# Calculate partial_firmclaims_partial_price based on candidate theta
		average_claims = exp.(X_CLMS * theta);
		ac_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * average_claims;
		demand_partial_ac_partial_price = theta[findall(claims_moments_covariates .== "log_risk_score")[1]] .*
			(transpose(rs_partial_matrix .* own_plan_year_matrix) * (average_claims .*
				plan_demands ./ plan_risk_scores));
		partial_firmclaims_partial_price_theta = (1 ./ partial_plan_demand_partial_price) .* 
			(ac_partial_ownplans_partial_price .+ demand_partial_ac_partial_price);
		
		# Calculate FOC values 
		foc_values = (transpose(foc_moment_indicator) * 
				((partial_firmclaims_partial_price_focs .- partial_firmclaims_partial_price_theta) .* plan_demands)) ./ 
			(transpose(foc_moment_indicator) * plan_demands);
		
		if m > 0
			foc_values = foc_values[m];
		end
		
		return(foc_values)
	end
	
	function compute_avg_claims_moment_value(theta_CLMS,m)
		avg_claims_moment_values = transpose(X_claims_data .* (ac_weights ./ sum(ac_weights))) * 
			(exp.(y_claims_data) .- exp.(X_claims_data * theta_CLMS));
		return(avg_claims_moment_values[m])
	end

# Calculate the key partial derivatives	
		
	function calculate_key_partials(theta_RS,pmt,hdata,X_RS,lambda)
		
		# Setup variables
		base_premium = plans[!,:Premium];
		probabilities = hdata[:,findall(hdata_fields .== "real_prob")[1]];
		ex_probs = hdata[:,findall(hdata_fields .== "real_ex_prob")[1]];
		alpha = hdata[:,findall(hdata_fields .== "alpha")[1]];
		rating_factor = hdata[:,findall(hdata_fields .== "rating_factor")[1]];
		slc_plans = hdata[:,findall(hdata_fields .== "real_slc_plan")[1]];
		slc_indices = findall(hdata[:,findall(hdata_fields .== "real_slc_plan")[1]] .== 1);
		
		# Calculate weights
		probability_weights = probabilities .* hdata[:,findall(hdata_fields .== "weight")[1]];
		pricing_weights = probabilities .* hdata[:,findall(hdata_fields .== "age_factor")[1]];
		plan_demands = pmt[:,findall(pmt_fields .== "real_demand")[1]];
		plan_pricing_weights = pmt[:,findall(pmt_fields .== "plan_pricing_weight")[1]];
		
		# Update risk score variables
		plan_risk_scores = exp.(X_RS * theta_RS);
		theta_RS_share_variables = theta_RS[theta_RS_share_indices];
		household_risk_scores = plan_risk_scores[household_pmt_indices];
			
		# Calculate cost factors
		rs_cost_factor = plans_pmt[!,:IDF] .* plans_pmt[!,:GCF] .* plan_risk_scores;
		ra_cost_factor = plans_pmt[!,:AV] .* plans_pmt[!,:IDF] .* plans_pmt[!,:GCF] .* pmt[:,findall(pmt_fields .== "average_age_factor")[1]];
		household_ra_cost_factor = (hdata[:,findall(hdata_fields .== "age_factor")[1]] ./ hdata[:,findall(hdata_fields .== "members")[1]]) .*
			hdata[:,findall(hdata_fields .== "ra_factor")[1]] .* hdata[:,findall(hdata_fields .== "GCF")[1]];
		
		# Initialize household-level variables that depend on household jacobian
		jacobian = zeros(length(plan_pmt_names),length(plan_pmt_names));
		rs_partial_matrix = zeros(length(plan_pmt_names),length(plan_pmt_names));
		partials_to_return = zeros(length(plan_pmt_names),length(partial_names));
		price_partial_firmprob_partial_price_matrix = zeros(length(insurers),length(plan_pmt_names));
		
		partial_rs_demographic_own_demand_partial_price = zeros(length(plan_pmt_names));
		partial_rs_demographic_all_demand_partial_price = zeros(length(plan_pmt_names));
		
		partial_ra_demographic_own_demand_partial_price = zeros(length(plan_pmt_names));
		partial_ra_demographic_all_demand_partial_price = zeros(length(plan_pmt_names));
		
		# Base vectors to form variables that depend on household-level jacobian
		cross_vector1 = (lambda-1)/lambda * ex_probs .- probabilities;
		cross_vector2 = probabilities .* alpha;
		own_vector = (1/lambda .+ (lambda-1)/lambda * ex_probs .- probabilities) .* probabilities .* alpha;
		household_premiums = rating_factor .* base_premium[household_pt_indices];
		
		# NOTE: You can make this loop more modular, but it slows down the code
		
		for m in 1:size(unique_choice_sets)[1]
			
			# Indices
		
				# Get data indices
				choice_set_indices = findall(data[uninsured_plans .== 0,:choice_set_id] .== m);
				subsidized_indices = intersect(choice_set_indices,findall(households[household_numbers,:subsidized_members] .> 0));
				unsubsidized_indices = intersect(choice_set_indices,findall(households[household_numbers,:subsidized_members] .== 0));
					
				# Get household indices	
				household_indices = unique(household_numbers[choice_set_indices]);
				subsidized_household_indices = unique(household_numbers[subsidized_indices]);
				unsubsidized_household_indices = unique(household_numbers[unsubsidized_indices]); 
				
				# Get plan indices
				jacobian_indices = unique(household_pmt_indices[choice_set_indices]);
				if length(subsidized_indices) > 0
					choice_slc_jacobian_index = findall(slc_plans[subsidized_indices] .> 0)[1];
					choice_nonslc_jacobian_indices = setdiff(1:length(jacobian_indices),choice_slc_jacobian_index);
				end
				firm_indices = findall(vec(sum(dims=2,firm_plan_matrix[:,jacobian_indices])) .> 0);	
					
				num_households_subs = length(subsidized_household_indices);
				num_households_unsubs = length(unsubsidized_household_indices);
				num_plans = length(jacobian_indices);
				
			# Matrices we will reuse many times to calculate all jacobian variants
			if num_households_unsubs > 0
				cross_matrix1_unsubs = reshape(cross_vector1[unsubsidized_indices],num_plans,num_households_unsubs);
				cross_matrix2_unsubs = reshape(cross_vector2[unsubsidized_indices],num_plans,num_households_unsubs);
				own_matrix_unsubs = reshape(own_vector[unsubsidized_indices],num_plans,num_households_unsubs);
				household_premiums_unsubs = reshape(household_premiums[unsubsidized_indices],num_plans,num_households_unsubs);
			end
			if num_households_subs > 0
				cross_matrix1_subs = reshape(cross_vector1[subsidized_indices],num_plans,num_households_subs);
				cross_matrix2_subs = reshape(cross_vector2[subsidized_indices],num_plans,num_households_subs);
				own_matrix_subs = reshape(own_vector[subsidized_indices],num_plans,num_households_subs);
				household_premiums_subs = reshape(household_premiums[subsidized_indices],num_plans,num_households_subs);
			end	
				
			# Get price jacobian
			
				if num_households_unsubs > 0
					multiplier_unsubs = reshape(rating_factor[unsubsidized_indices] .* household_weights[unsubsidized_indices],num_plans,num_households_unsubs);
					unsubsidized_portion = cross_matrix1_unsubs * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs);
					jacobian[jacobian_indices,jacobian_indices] += unsubsidized_portion;
				end
				
				if num_households_subs > 0
					multiplier_subs = reshape(rating_factor[subsidized_indices] .* household_weights[subsidized_indices],num_plans,num_households_subs);
					subsidized_portion = cross_matrix1_subs * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs);
					
					# Apply adjustment for endogenous subsidy (only affects SLC column of subsidized portion)
							# want to calculate partial q_l/partial p_SLC for all l (including SLC)
							# negative sum of all partial q_l/partial p_l', except l' = SLC 
					
					multiplier_slc = reshape(repeat(rating_factor[intersect(subsidized_indices,slc_indices)],
						inner=length(jacobian_indices)) .* household_weights[subsidized_indices],num_plans,num_households_subs);
					slc_portion = cross_matrix1_subs * transpose(cross_matrix2_subs .* multiplier_slc);
					slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc);
					subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices]));

					jacobian[jacobian_indices,jacobian_indices] += subsidized_portion;
				end
				
			# Populate the risk score partial matrix (input for the claim slope matrix)
				# Note this is a portion; we'll subtract the second portion later
			
				if num_households_unsubs > 0
					rs_multiplier_unsubs = reshape((rs_cost_factor[jacobian_indices] ./ plan_demands[jacobian_indices]) * 
						transpose(hh_share_matrix[unsubsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_unsubs);
					unsubsidized_portion = (cross_matrix1_unsubs .* rs_multiplier_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* rs_multiplier_unsubs);
					rs_partial_matrix[jacobian_indices,jacobian_indices] += unsubsidized_portion;
				end
				
				if num_households_subs > 0
					rs_multiplier_subs = reshape((rs_cost_factor[jacobian_indices] ./ plan_demands[jacobian_indices]) * 
						transpose(hh_share_matrix[subsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_subs);
					subsidized_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* rs_multiplier_subs);
						
					slc_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
					slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* rs_multiplier_subs);
					subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices])); 	

					rs_partial_matrix[jacobian_indices,jacobian_indices] += subsidized_portion;						
				end
				
				
			# Get partials that go into marginal RA transfer

				if num_households_unsubs > 0
					rs_multiplier_unsubs = reshape(rs_cost_factor[jacobian_indices] * 
						transpose(hh_share_matrix[unsubsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_unsubs);
					unsubsidized_portion = (cross_matrix1_unsubs .* rs_multiplier_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* rs_multiplier_unsubs);
				
					partial_rs_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_rs_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion));
					
					ra_multiplier_unsubs = reshape(household_ra_cost_factor[unsubsidized_indices],num_plans,num_households_unsubs);;
					unsubsidized_portion = (cross_matrix1_unsubs .* ra_multiplier_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* ra_multiplier_unsubs);
					
					partial_ra_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_ra_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,unsubsidized_portion));
				
				end

				if num_households_subs > 0
					rs_multiplier_subs = reshape(rs_cost_factor[jacobian_indices] * 
						transpose(hh_share_matrix[subsidized_household_indices,:] * theta_RS_share_variables),num_plans,num_households_subs);
					subsidized_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* rs_multiplier_subs);
					
					slc_portion = (cross_matrix1_subs .* rs_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
					slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* rs_multiplier_subs);
					subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices]));
					
					partial_rs_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_rs_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion));

					ra_multiplier_subs = reshape(household_ra_cost_factor[subsidized_indices],num_plans,num_households_subs);;
					subsidized_portion = (cross_matrix1_subs .* ra_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* ra_multiplier_subs);

					slc_portion = (cross_matrix1_subs .* ra_multiplier_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
					slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* ra_multiplier_subs);
					subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices]));

					partial_ra_demographic_own_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion .* own_plan_year_matrix[jacobian_indices,jacobian_indices]));
					partial_ra_demographic_all_demand_partial_price[jacobian_indices] += vec(sum(dims=1,subsidized_portion));
				end
			
			# Get partial_firmprob_partial_price (do not multiply by household's weight)
			
				if num_households_unsubs > 0
					multiplier_unsubs = reshape(rating_factor[unsubsidized_indices],num_plans,num_households_unsubs);
					unsubsidized_portion = (cross_matrix1_unsubs .* household_premiums_unsubs) * transpose(cross_matrix2_unsubs .* multiplier_unsubs);
					unsubsidized_portion[diagind(unsubsidized_portion)] = sum(dims=2,own_matrix_unsubs .* multiplier_unsubs .* household_premiums_unsubs);
					price_partial_firmprob_partial_price_matrix[firm_indices,jacobian_indices] += 
						firm_plan_matrix[firm_indices,jacobian_indices] * unsubsidized_portion;
				end

				if num_households_subs > 0
					multiplier_subs = reshape(rating_factor[subsidized_indices],num_plans,num_households_subs);
					subsidized_portion = (cross_matrix1_subs .* household_premiums_subs) * transpose(cross_matrix2_subs .* multiplier_subs);
					subsidized_portion[diagind(subsidized_portion)] = sum(dims=2,own_matrix_subs .* multiplier_subs .* household_premiums_subs);
					
					choice_slc_indices = intersect(subsidized_indices,slc_indices);
					multiplier_slc = reshape(repeat(rating_factor[choice_slc_indices],inner=length(jacobian_indices)),num_plans,num_households_subs);
					slc_portion = (cross_matrix1_subs .* household_premiums_subs) * transpose(cross_matrix2_subs .* multiplier_slc);
					slc_portion[diagind(slc_portion)] = sum(dims=2,own_matrix_subs .* multiplier_slc .* household_premiums_subs);
					subsidized_portion[:,choice_slc_jacobian_index] = -vec(sum(dims=2,slc_portion[:,choice_nonslc_jacobian_indices])); 
					
					price_partial_firmprob_partial_price_matrix[firm_indices,jacobian_indices] += 
						firm_plan_matrix[firm_indices,jacobian_indices] * subsidized_portion;
				end
		end
		
		# Calculate marginal revenue
			
		partial_firmrevenue_partial_price_matrix = price_partial_firmprob_partial_price_matrix .+
			firm_plan_matrix .* transpose(plan_pricing_weights);
		own_firmrevenue_indices = findall(transpose(partial_firmrevenue_partial_price_matrix .* firm_plan_matrix) .!= 0);			
		partials_to_return[:,findall(partial_names .== "partial_firmrevenue_partial_price")[1]] = 
			transpose(partial_firmrevenue_partial_price_matrix)[own_firmrevenue_indices];
		partial_marketrevenue_partial_price = vec(sum(dims=1,partial_firmrevenue_partial_price_matrix));
				
		# Calculate marginal transfer
		
		annual_total_premiums = transpose(plan_year_matrix) * pmt[:,findall(pmt_fields .== "total_premiums")[1]];
		annual_demand = transpose(plan_year_matrix) * plan_demands;
		annual_ra_weights = transpose(plan_year_matrix) * (plan_demands .* ra_cost_factor);
		annual_rs_weights = transpose(plan_year_matrix) * (plan_demands .* rs_cost_factor);
			
		insurer_rs_weights = transpose(own_plan_year_matrix) * (plan_demands .* rs_cost_factor);
		insurer_ra_weights = transpose(own_plan_year_matrix) * (plan_demands .* ra_cost_factor);	
		insurer_rs_shares = insurer_rs_weights ./ annual_rs_weights;
		insurer_ra_shares = insurer_ra_weights ./ annual_ra_weights;
		
		rs_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * rs_cost_factor;
		rs_partial_allplans_partial_price = transpose(jacobian) * rs_cost_factor;
		ra_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * ra_cost_factor;
		ra_partial_allplans_partial_price = transpose(jacobian) * ra_cost_factor;
		
		partial_rs_dem_share_own_plan_demand_partial_price = transpose(jacobian .* own_plan_year_matrix) * 
			(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables));
		partial_rs_dem_share_all_plan_demand_partial_price = transpose(jacobian) * 
			(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables));		
				partial_ra_dem_share_own_plan_demand_partial_price = (transpose(jacobian .* own_plan_year_matrix) * ra_cost_factor);
		partial_ra_dem_share_all_plan_demand_partial_price = (transpose(jacobian) * ra_cost_factor);	
				own_plan_firmshare_partial_rs_partial_price =
			partial_rs_demographic_own_demand_partial_price .- partial_rs_dem_share_own_plan_demand_partial_price;
		own_plan_firmshare_partial_ra_partial_price =
			partial_ra_demographic_own_demand_partial_price .- partial_ra_dem_share_own_plan_demand_partial_price;
		
		all_plan_firmshare_partial_rs_partial_price =			
			partial_rs_demographic_all_demand_partial_price .- partial_rs_dem_share_all_plan_demand_partial_price;
		all_plan_firmshare_partial_ra_partial_price =			
			partial_ra_demographic_all_demand_partial_price .- partial_ra_dem_share_all_plan_demand_partial_price;
			
		partial_firmriskshare_partial_price = (1 ./ annual_rs_weights) .*
			((own_plan_firmshare_partial_rs_partial_price .+ rs_partial_ownplans_partial_price) .-
				insurer_rs_shares .* (all_plan_firmshare_partial_rs_partial_price .+ rs_partial_allplans_partial_price));	
		partial_firmutilizationshare_partial_price = (1 ./ annual_ra_weights) .*
			((own_plan_firmshare_partial_ra_partial_price .+ ra_partial_ownplans_partial_price) .-
				insurer_ra_shares .* (all_plan_firmshare_partial_ra_partial_price .+ ra_partial_allplans_partial_price));	
		
		partials_to_return[:,findall(partial_names .== "partial_firmtransfer_partial_price")[1]] = nu_for_riskadj .*
			(annual_total_premiums .* (partial_firmriskshare_partial_price .- partial_firmutilizationshare_partial_price) .+
			partial_marketrevenue_partial_price .* (insurer_rs_shares .- insurer_ra_shares));
			
		# Calculate marginal admin
		partials_to_return[:,findall(partial_names .== "partial_firmadmin_partial_price")[1]] =
			transpose(jacobian .* own_plan_year_matrix) * pmt[:,findall(pmt_fields .== "variable_admin")[1]];
			
		# Calculate second part of risk score partial matrix
		rs_partial_matrix -= jacobian .*
			(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables) ./ plan_demands);
		
		# Calculate partial firm demand partial plan demand 
		partials_to_return[:,findall(partial_names .== "partial_plan_demand_partial_price")[1]] = diag(jacobian);
			
		return(cat(dims=2,jacobian,rs_partial_matrix,partials_to_return))# Calculate marginal revenue
			
		partial_firmrevenue_partial_price_matrix = price_partial_firmprob_partial_price_matrix .+
			firm_plan_matrix .* transpose(plan_pricing_weights);
		own_firmrevenue_indices = findall(transpose(partial_firmrevenue_partial_price_matrix .* firm_plan_matrix) .!= 0);			
		partials_to_return[:,findall(partial_names .== "partial_firmrevenue_partial_price")[1]] = 
			transpose(partial_firmrevenue_partial_price_matrix)[own_firmrevenue_indices];
		partial_marketrevenue_partial_price = vec(sum(dims=1,partial_firmrevenue_partial_price_matrix));
				
		# Calculate marginal transfer
		
			annual_total_premiums = transpose(plan_year_matrix) * pmt[:,findall(pmt_fields .== "total_premiums")[1]];
			annual_demand = transpose(plan_year_matrix) * plan_demands;
			annual_ra_weights = transpose(plan_year_matrix) * (plan_demands .* ra_cost_factor);
			annual_rs_weights = transpose(plan_year_matrix) * (plan_demands .* rs_cost_factor);
			
			insurer_rs_weights = transpose(own_plan_year_matrix) * (plan_demands .* rs_cost_factor);
			insurer_ra_weights = transpose(own_plan_year_matrix) * (plan_demands .* ra_cost_factor);	
			insurer_rs_shares = insurer_rs_weights ./ annual_rs_weights;
			insurer_ra_shares = insurer_ra_weights ./ annual_ra_weights;
			
			rs_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * rs_cost_factor;
			rs_partial_allplans_partial_price = transpose(jacobian) * rs_cost_factor;
			ra_partial_ownplans_partial_price = transpose(jacobian .* own_plan_year_matrix) * ra_cost_factor;
			ra_partial_allplans_partial_price = transpose(jacobian) * ra_cost_factor;
			
			partial_rs_dem_share_own_plan_demand_partial_price = transpose(jacobian .* own_plan_year_matrix) * 
				(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables));
			partial_rs_dem_share_all_plan_demand_partial_price = transpose(jacobian) * 
				(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables));		
	
			partial_ra_dem_share_own_plan_demand_partial_price = (transpose(jacobian .* own_plan_year_matrix) * ra_cost_factor);
			partial_ra_dem_share_all_plan_demand_partial_price = (transpose(jacobian) * ra_cost_factor);	
	
			own_plan_firmshare_partial_rs_partial_price =
				partial_rs_demographic_own_demand_partial_price .- partial_rs_dem_share_own_plan_demand_partial_price;
			own_plan_firmshare_partial_ra_partial_price =
				partial_ra_demographic_own_demand_partial_price .- partial_ra_dem_share_own_plan_demand_partial_price;
			
			all_plan_firmshare_partial_rs_partial_price =			
				partial_rs_demographic_all_demand_partial_price .- partial_rs_dem_share_all_plan_demand_partial_price;
			all_plan_firmshare_partial_ra_partial_price =			
				partial_ra_demographic_all_demand_partial_price .- partial_ra_dem_share_all_plan_demand_partial_price;
			
			partial_firmriskshare_partial_price = (1 ./ annual_rs_weights) .*
				((own_plan_firmshare_partial_rs_partial_price .+ rs_partial_ownplans_partial_price) .-
					insurer_rs_shares .* (all_plan_firmshare_partial_rs_partial_price .+ rs_partial_allplans_partial_price));	
			partial_firmutilizationshare_partial_price = (1 ./ annual_ra_weights) .*
				((own_plan_firmshare_partial_ra_partial_price .+ ra_partial_ownplans_partial_price) .-
					insurer_ra_shares .* (all_plan_firmshare_partial_ra_partial_price .+ ra_partial_allplans_partial_price));	
			
			partials_to_return[:,findall(partial_names .== "partial_firmtransfer_partial_price")[1]] = nu_for_riskadj .*
				(annual_total_premiums .* (partial_firmriskshare_partial_price .- partial_firmutilizationshare_partial_price) .+
				partial_marketrevenue_partial_price .* (insurer_rs_shares .- insurer_ra_shares));
			
		# Calculate marginal admin
		partials_to_return[:,findall(partial_names .== "partial_firmadmin_partial_price")[1]] =
			transpose(jacobian .* own_plan_year_matrix) * pmt[:,findall(pmt_fields .== "variable_admin")[1]];
			
		# Calculate second part of risk score partial matrix
		rs_partial_matrix -= jacobian .*
			(rs_cost_factor .* (X_RS[:,theta_RS_share_indices] * theta_RS_share_variables) ./ plan_demands);
		
		# Calculate partial firm demand partial plan demand 
		partials_to_return[:,findall(partial_names .== "partial_plan_demand_partial_price")[1]] = diag(jacobian);
			
		return(cat(dims=2,jacobian,rs_partial_matrix,partials_to_return))
	end			
		
	
# Determine SLC plan for each household

	function determine_SLC_plan(silver_premiums)
	
		if length(silver_premiums) == 1
			slc_silver = 1;
		elseif length(silver_premiums) > 1 
			cheapest_silver = argmin(silver_premiums);
			slc_silver = findall(silver_premiums .== minimum(silver_premiums[setdiff(1:length(silver_premiums),cheapest_silver)]))[1];
		else
			println("Error: No silver plans!")
		end
	
		return(slc_silver)
	end

# Initialize hdata object (household data)	
	
	function initialize_household_data(X,hdata_fields,data,plans,household_numbers,min_beta)
		
		if nested
			lambda = last(min_beta);
			min_beta = min_beta[1:(length(min_beta)-1)];
		end
		
		hdata = zeros(Float64,length(household_numbers),length(hdata_fields)) ;
		
		hdata[:,findall(hdata_fields .== "penalty")[1]] = convert(Array,data[uninsured_plans .== 0,:penalty]);
		hdata[:,findall(hdata_fields .== "weight")[1]] = convert(Array,data[uninsured_plans .== 0,:weight]);
		hdata[:,findall(hdata_fields .== "choice")[1]] = convert(Array,data[uninsured_plans .== 0,:choice]);
		hdata[:,findall(hdata_fields .== "pricing_factor")[1]] = convert(Array,data[uninsured_plans .== 0,:pricing_factor]);
		hdata[:,findall(hdata_fields .== "rating_factor")[1]] = convert(Array,data[uninsured_plans .== 0,:age_geog_factor]);
		hdata[:,findall(hdata_fields .== "age_factor")[1]] = convert(Array,data[uninsured_plans .== 0,:age_factor]);
		hdata[:,findall(hdata_fields .== "ra_factor")[1]] = convert(Array,data[uninsured_plans .== 0,:ra_factor]);
		hdata[:,findall(hdata_fields .== "GCF")[1]] = convert(Array,data[uninsured_plans .== 0,:GCF]);
		hdata[:,findall(hdata_fields .== "alpha")[1]] = calculate_alpha(min_beta,household_numbers,"All");
		hdata[:,findall(hdata_fields .== "members")[1]] .= households[household_numbers,:enrollees];				
				
		# Household unsubsidized premiums
		for j in plan_year_names
			plan_indices = findall(data_pt_names .== j);			
			hdata[plan_indices,findall(hdata_fields .== "real_unsubs_premium")[1]] = 
				plans[findall(plans[!,:pt_name] .== j)[1],:Premium] .* 
				hdata[plan_indices,findall(hdata_fields .== "rating_factor")[1]];
		end
		
		
		# SLC Premium
		
			# Create data frame to hold objects we will group by household
			silver_indices = sort(cat(dims=1,findall(data[uninsured_plans .== 0,:AV] .== .7),findall(data[uninsured_plans .== 0,:AV] .== .73),
				findall(data[uninsured_plans .== 0,:AV] .== .87),findall(data[uninsured_plans .== 0,:AV] .== .94)));
			silverdf = DataFrame(household_number = household_numbers[silver_indices],
				premium = hdata[silver_indices,findall(hdata_fields .== "real_unsubs_premium")[1]],silver_indices = silver_indices);
			
			# Determine SLC plans by household
			num_household_silver_plans = combine(groupby(silverdf,:household_number), :household_number => length)[!,:household_number_length];
			pop!(num_household_silver_plans)
			num_household_silver_plans = cat(dims=1,0,num_household_silver_plans);
			starting_indices = cumsum(num_household_silver_plans); # starting position for each household
			household_slc_indices = combine(groupby(silverdf,:household_number), :premium => determine_SLC_plan)[!,:premium_determine_SLC_plan];
			slc_indices = silverdf[starting_indices .+ household_slc_indices,:silver_indices]; # SLC indices in hdata
			hdata[slc_indices,findall(hdata_fields .== "real_slc_plan")[1]] = (households[!,:subsidized_members] .> 0); # SLC index gets 1 if house is subsidized, 0 otherwise
		
		# Choice probabilities
		for n in 1:ns
			X_demand = copy(X);
			X_demand[:,premium_indices] = X_demand[:,premium_indices] .* X_demand[:,findall(all_covariates .== "Premium")[1]];;
			X_demand[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
			
			expdf = DataFrame(household_number = household_numbers,expvarall = exp.(X_demand[uninsured_plans .== 0,:] * min_beta));
			expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			expsum = max.(expsum,0.0000000000001); # if expsum is computationally 0, then the uninsured probability will be exactly 1, which causes ex_probs to be infinite
			
			
			if inattention
				entrants = findall(households[!,:previous_plan_number] .== "NA");
				search_probabilities = 1 ./ (1 .+ exp.(Xsearch * beta_search));
				search_probabilities[entrants] .= 1; 
				
				uninsured_probabilities = 1  ./ (1 .+ expsum);
				insured_probabilities = search_probabilities[household_numbers] .* expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers]) .+ 
					(1 .- search_probabilities[household_numbers]) .* data[uninsured_plans .== 0,:previous_choice];
				hdata[:,findall(hdata_fields .== "real_prob")[1]] += insured_probabilities ./ ns;
				hdata[:,findall(hdata_fields .== "real_ex_prob")[1]] += (insured_probabilities ./ (1 .- uninsured_probabilities[household_numbers])) ./ ns;
			else 
				uninsured_probabilities = 1  ./ (1 .+ expsum);
				insured_probabilities = expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers]);
				hdata[:,findall(hdata_fields .== "real_prob")[1]] += insured_probabilities ./ ns;
				hdata[:,findall(hdata_fields .== "real_ex_prob")[1]] += (insured_probabilities ./ (1 .- uninsured_probabilities[household_numbers])) ./ ns;
				
				# Instrumented prob
				X_demand[:,findall(all_covariates .== "Premium")[1]] .-= X_demand[:,findall(all_covariates .== "residual")[1]];
				expdf = DataFrame(household_number = household_numbers,expvarall = exp.(X_demand[uninsured_plans .== 0,:] * min_beta));
				expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
				hdata[:,findall(hdata_fields .== "real_instrumented_prob")[1]] += (expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers])) ./ ns;
			end
		
			
		end
			
		return(hdata)
	end

# Rescale weights

	function rescale_weights(moments_data,weights,yearmin,yearmax)
		
		# This function weights within insurer, but assigns each insurer the same total weight
		
		rescaled_weights = zeros(length(weights));
		year_indices = intersect(findall(moments_data[!,:year] .>= yearmin),findall(moments_data[!,:year] .<= yearmax));
		small_insurer_plans = 1 .- (moments_data[year_indices,:Anthem] .+ moments_data[year_indices,:Blue_Shield] .+ 
			moments_data[year_indices,:Health_Net] .+ moments_data[year_indices,:Kaiser]);
		
		for f in insurers
			for t in yearmin:yearmax
				if f == "Small_Insurer"
					insurer_year_indices = intersect(findall(small_insurer_plans .== 1),
						findall(moments_data[year_indices,:year] .== t));
				else
					insurer_year_indices = intersect(findall(moments_data[year_indices,findall(names(moments_data) .== f)[1]] .== 1),
						findall(moments_data[year_indices,:year] .== t));
				end
				rescaled_weights[insurer_year_indices] = (weights[insurer_year_indices] ./ sum(weights[insurer_year_indices]))/length(insurers);
			end
		end
		
		return(rescaled_weights)
	end

# Create household demographic share matrix

	function update_dem_share_matrix(households)
	
		hh_share_matrix = zeros(size(households)[1],length(risk_score_share_variables));
		
		if all_rs_variables 
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_0to34")[1]] = households[!,:perc_0to17] .+ households[!,:perc_18to25] .+ households[!,:perc_26to34];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_35to54")[1]] = households[!,:perc_35to44] .+ households[!,:perc_45to54];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_250to400")[1]] = ((households[!,:FPL] .> 2.50) .& (households[!,:FPL] .<= 4)) .* 1;
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_gt400")[1]] = (households[!,:FPL] .> 4) .* 1;
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_male")[1]] = households[!,:perc_male];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_family")[1]] = (households[!,:household_size] .> 1) .* 1;
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_Asian")[1]] = households[!,:perc_asian];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_Black")[1]] = households[!,:perc_black];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_Hispanic")[1]] = households[!,:perc_hispanic];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_Other_Race")[1]] = households[!,:perc_other];
		elseif int_rs_variables
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_0to34")[1]] = households[!,:perc_0to17] .+ households[!,:perc_18to25] .+ households[!,:perc_26to34];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_male")[1]] = households[!,:perc_male];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_family")[1]] = (households[!,:household_size] .> 1) .* 1;
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_minority")[1]] = 
				households[!,:perc_asian] .+ households[!,:perc_black] .+ households[!,:perc_hispanic] .+ households[!,:perc_other];
		elseif hisp_flag
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_0to34")[1]] = households[!,:perc_0to17] .+ households[!,:perc_18to25] .+ households[!,:perc_26to34];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_male")[1]] = households[!,:perc_male];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_family")[1]] = (households[!,:household_size] .> 1) .* 1;
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_Hispanic")[1]] = households[!,:perc_hispanic];
		elseif agegender_flag
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_0to34")[1]] = households[!,:perc_0to17] .+ households[!,:perc_18to25] .+ households[!,:perc_26to34];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_male")[1]] = households[!,:perc_male];
		elseif old_flag
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_18to54")[1]] = 
				households[!,:perc_18to25] .+ households[!,:perc_26to34] .+ households[!,:perc_35to44] .+ households[!,:perc_45to54];
			hh_share_matrix[:,findall(risk_score_share_variables .== "share_Hispanic")[1]] = households[!,:perc_hispanic];
		end
		
		return(hh_share_matrix)
	end

# Initialize pmt object (plan data at the plan-market-year level)		
	
	function initialize_pmt_data(pmt_fields,plans_pmt,hdata)

		pmt = zeros(size(plans_pmt)[1],length(pmt_fields));
		
		weight_object = DataFrame(household_number = household_numbers, plan_name = data_pmt_names,
							probability_weights = hdata[:,findall(hdata_fields .== "real_prob")[1]] .* hdata[:,findall(hdata_fields .== "weight")[1]],
							pricing_weights = hdata[:,findall(hdata_fields .== "age_factor")[1]] .* hdata[:,findall(hdata_fields .== "real_prob")[1]],
							plan_pricing_weights = hdata[:,findall(hdata_fields .== "rating_factor")[1]] .* hdata[:,findall(hdata_fields .== "real_prob")[1]],
							age_weights = (hdata[:,findall(hdata_fields .== "age_factor")[1]] ./ hdata[:,findall(hdata_fields .== "members")[1]]) .* 
								hdata[:,findall(hdata_fields .== "real_prob")[1]] .* hdata[:,findall(hdata_fields .== "weight")[1]]);
		
		
		pmt[:,findall(pmt_fields .== "real_demand")[1]] = 
			sort!(combine(groupby(weight_object,:plan_name), :probability_weights => sum),[:plan_name])[!,:probability_weights_sum];
		pmt[:,findall(pmt_fields .== "plan_pricing_weight")[1]] = 
			sort!(combine(groupby(weight_object,:plan_name), :pricing_weights => sum),[:plan_name])[!,:pricing_weights_sum];
		pmt[:,findall(pmt_fields .== "total_premiums")[1]] = plans[year_plan_indices,:Premium] .*
			sort!(combine(groupby(weight_object,:plan_name), :plan_pricing_weights => sum),[:plan_name])[!,:plan_pricing_weights_sum];
		pmt[:,findall(pmt_fields .== "average_age_factor")[1]] = 
			sort!(combine(groupby(weight_object,:plan_name), :age_weights => sum),[:plan_name])[!,:age_weights_sum] ./
			pmt[:,findall(pmt_fields .== "real_demand")[1]];		
			
		pmt[:,findall(pmt_fields .== "reins_factor")[1]] = plans[year_plan_indices,:reins_factor];
		pmt[:,findall(pmt_fields .== "variable_admin")[1]] = plans[year_plan_indices,:variable_admin];
		pmt[:,findall(pmt_fields .== "fixed_admin")[1]] = plans[year_plan_indices,:fixed_admin];
		pmt[:,findall(pmt_fields .== "base_premium")[1]] = plans[year_plan_indices,:Premium];
		
		# RA factors (from CMS RA transfer paper)
		pmt[findall(plans_pmt[!,:AV] .== 0.6),findall(pmt_fields .== "ra_factor")[1]] .= 0.6;
		pmt[findall(plans_pmt[!,:AV] .== 0.7),findall(pmt_fields .== "ra_factor")[1]] .= 0.7 * 1.03;
		pmt[findall(plans_pmt[!,:AV] .== 0.8),findall(pmt_fields .== "ra_factor")[1]] .= 0.8 * 1.08;
		pmt[findall(plans_pmt[!,:AV] .== 0.9),findall(pmt_fields .== "ra_factor")[1]] .= 0.9 * 1.15;
		
		return(pmt)
	end
	
# Calculate Price Coefficient

	# NOTE: The previous choice interaction will need to be removed for 2014 only

	function calculate_alpha(betap,household_numbers,skip) 
	
		# Need to get the premium coefficient alpha for each household
		
		alpha = ones(length(household_numbers)) .* betap[findall(all_covariates .== "Premium")[1]];
		
		# Income
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_250to400")[1]] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) ; 	
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_gt400")[1]] .* (households[household_numbers,:FPL] .> 4) ;
			
		# Age
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_0to34")[1]] .* (households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_35to54")[1]] .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
		
		if add_complex_interactions
			alpha = alpha .+ betap[findall(all_covariates .== "Premium_250to400_0to34")[1]] .*
				((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
				(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
			alpha = alpha .+ betap[findall(all_covariates .== "Premium_gt400_0to34")[1]] .*
				((households[household_numbers,:FPL] .> 4)) .*
				(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
			alpha = alpha .+ betap[findall(all_covariates .== "Premium_250to400_35to54")[1]] .*
				((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
				(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
			alpha = alpha .+ betap[findall(all_covariates .== "Premium_gt400_35to54")[1]] .*
				((households[household_numbers,:FPL] .> 4)) .*
				(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
		end
		
		# Gender
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_male")[1]] .* (households[household_numbers,:perc_male]);
			
		# Family
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_family")[1]] .* (households[household_numbers,:household_size] .> 1);
		
		# Race
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_Black")[1]] .* (households[household_numbers,:perc_black]);
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_Hispanic")[1]] .* (households[household_numbers,:perc_hispanic]);	
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_Asian")[1]] .* (households[household_numbers,:perc_asian]);	
		alpha = alpha .+ betap[findall(all_covariates .== "Premium_Other")[1]] .* (households[household_numbers,:perc_other]);			
		
		return(alpha)
	end

# Form Covariate Matrix/Array X (Demand-Side Variables)

	function form_X_matrix(data,households,plans,product_covariates);
	
			X = zeros(Float64,size(data)[1],length(all_covariates)) ;
	
			household_indices = convert(Array,data[!,:household_number]);	
			uninsured_plans = convert(Array,data[!,:uninsured_plan]);
	
			# Add Subsidized Premium and Penalty (my premium variable includes subsidy and penalty)
			X[:,1] = convert(Array,data[!,:premium]);
		
			# AV is now an individual-specific variable to deal with CSR plans
			av_ind = findall(all_covariates .== "AV")[1];
			X[:,av_ind] = convert(Array,data[!,:AV]);
		
			data_pt_full_names = string.(data[!,:plan_name],data[!,:year]);
		
			for plan in plans[!,:pt_name]
				#plan_indices = findall(data[!,:plan_name] .== plan);
				plan_indices = findall(data_pt_full_names .== plan);
				for covariate in setdiff(product_covariates,["Premium","AV","Anthem_HMO"])
					plans_index_name = findall(plans[!,:pt_name] .== plan)[1];
					plans_index = findall(names(plans) .== covariate)[1];
					array_index = findall(all_covariates .== covariate)[1];
					X[plan_indices,array_index] .= plans[plans_index_name,plans_index];
				end
			end
			
			#X[:,findall(all_covariates .== "Anthem_HMO")[1]] = 
			#	X[:,findall(all_covariates .== "Anthem")[1]] .*
			#	X[:,findall(all_covariates .== "HMO")[1]];
			
			# Calculate per-capita premium
			X[:,1] = X[:,1] ./ households[household_indices,:enrollees];
			
			# Previous plan choice
			if "previous_choice" in all_covariates
				X[:,findall(all_covariates .== "previous_choice")[1]] = convert(Array,data[!,:previous_choice]);
				#X[:,findall(all_covariates .== "Premium_previous_choice")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]];
				X[:,findall(all_covariates .== "previous_choice_250to400")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* 
					((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) ; 
				X[:,findall(all_covariates .== "previous_choice_gt400")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* 
					(households[household_indices,:FPL] .> 4) ; 
				X[:,findall(all_covariates .== "previous_choice_0to34")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] + households[household_indices,:perc_26to34]);
				X[:,findall(all_covariates .== "previous_choice_35to54")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:perc_35to44] + households[household_indices,:perc_45to54]);
				X[:,findall(all_covariates .== "previous_choice_male")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:perc_male]);
				X[:,findall(all_covariates .== "previous_choice_family")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:household_size] .> 1);
				X[:,findall(all_covariates .== "previous_choice_Asian")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:perc_asian]);
				X[:,findall(all_covariates .== "previous_choice_Black")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:perc_black]);
				X[:,findall(all_covariates .== "previous_choice_Hispanic")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:perc_hispanic]);
				X[:,findall(all_covariates .== "previous_choice_Other")[1]] = 
					X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:perc_other]);
			
				if add_complex_interactions
					X[:,findall(all_covariates .== "previous_choice_250to400_0to34")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* ((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) .*
						(households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]); 
					X[:,findall(all_covariates .== "previous_choice_gt400_0to34")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:FPL] .> 4) .*
						(households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]);  
					X[:,findall(all_covariates .== "previous_choice_250to400_35to54")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* ((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) .*
						(households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]); 
					X[:,findall(all_covariates .== "previous_choice_gt400_35to54")[1]] = X[:,findall(all_covariates .== "previous_choice")[1]] .* (households[household_indices,:FPL] .> 4) .*
						(households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]);  
				end
				
				#X[:,findall(all_covariates .== "previous_choice_Anthem")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]] .*
				#	X[:,findall(all_covariates .== "Anthem")[1]];
				#X[:,findall(all_covariates .== "previous_choice_Blue_Shield")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]] .*
				#	X[:,findall(all_covariates .== "Blue_Shield")[1]];
				#X[:,findall(all_covariates .== "previous_choice_Kaiser")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]] .*
				#	X[:,findall(all_covariates .== "Kaiser")[1]];
				#X[:,findall(all_covariates .== "previous_choice_Health_Net")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]] .*
				#	X[:,findall(all_covariates .== "Health_Net")[1]];			
				#X[:,findall(all_covariates .== "previous_choice_HMO")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]] .*
				#	X[:,findall(all_covariates .== "HMO")[1]];
				#X[:,findall(all_covariates .== "previous_choice_AV")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]] .*
				#	X[:,findall(all_covariates .== "AV")[1]];
				#X[:,findall(all_covariates .== "previous_choice_silver")[1]] = 
				#	X[:,findall(all_covariates .== "previous_choice")[1]] .*
				#	X[:,findall(all_covariates .== "Silver")[1]];
				if network_flag 
					X[:,findall(all_covariates .== "previous_choice_breadth")[1]] = 
						X[:,findall(all_covariates .== "previous_choice")[1]] .*
						convert(Array,data[!,:network_breadth]);
					X[:,findall(all_covariates .== "previous_choice_inclus")[1]] = 
						X[:,findall(all_covariates .== "previous_choice")[1]] .*
						convert(Array,data[!,:network_inclus]);
						
					mean_breadth = sum(data[uninsured_plans .== 0,:network_breadth] .* data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight])/
						sum(data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight]);
					mean_inclus = sum(data[uninsured_plans .== 0,:network_inclus] .* data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight])/
						sum(data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight]);
					
					X[:,findall(all_covariates .== "breadth")[1]] = convert(Array,data[!,:network_breadth]) .- 
						mean_breadth * (1 .- uninsured_plans);
					X[:,findall(all_covariates .== "inclus")[1]] = convert(Array,data[!,:network_inclus]) .- 
						mean_inclus * (1 .- uninsured_plans);
					X[:,findall(all_covariates .== "pairinclus")[1]] = convert(Array,data[!,:network_pairinclus]);
				end
			end
			
			# Intercepts
			X[:,findall(all_covariates .== "intercept")] = 1 .- uninsured_plans;
			if !eq_robust_flag
				X[:,findall(all_covariates .== "intercept_year2015")] = (1 .- uninsured_plans) .* (households[household_indices,:year] .== 2015);
				X[:,findall(all_covariates .== "intercept_year2016")] = (1 .- uninsured_plans) .* (households[household_indices,:year] .== 2016);
			end
			X[:,findall(all_covariates .== "intercept_year2017")] = (1 .- uninsured_plans) .* (households[household_indices,:year] .== 2017);
			X[:,findall(all_covariates .== "intercept_year2018")] = (1 .- uninsured_plans) .* (households[household_indices,:year] .== 2018);
			X[:,findall(all_covariates .== "HMO_0to34")[1]] = X[:,findall(all_covariates .== "HMO")[1]] .* (households[household_indices,:perc_0to17] .+
				households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]);
			X[:,findall(all_covariates .== "HMO_35to54")[1]] = X[:,findall(all_covariates .== "HMO")[1]] .* 
				(households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]);
			
			#X[:,findall(all_covariates .== "intercept_male")] = (1 .- uninsured_plans) .*
			#	households[household_indices,:perc_male];
			#X[:,findall(all_covariates .== "intercept_family")] = (1 .- uninsured_plans) .*
			#	(households[household_indices,:household_size] .> 1);
			#X[:,findall(all_covariates .== "intercept_0to17")] = (1 .- uninsured_plans) .*
			#	(households[household_indices,:perc_0to17]);
			#X[:,findall(all_covariates .== "intercept_18to34")] = (1 .- uninsured_plans) .* 
			#	(households[household_indices,:perc_18to25] + households[household_indices,:perc_26to34]);	
			#X[:,findall(all_covariates .== "intercept_35to54")] = (1 .- uninsured_plans) .* 
			#	(households[household_indices,:perc_35to44] + households[household_indices,:perc_45to54]);
			#X[:,findall(all_covariates .== "intercept_250to400")] = (1 .- uninsured_plans) .* 
			#	((households[household_indices,:FPL] .> 2.5) .& (households[household_indices,:FPL] .<= 4)) ;						
			#X[:,findall(all_covariates .== "intercept_gt400")] = (1 .- uninsured_plans) .* 
			#	(households[household_indices,:FPL] .> 4) ; 
			#X[:,findall(all_covariates .== "intercept_asian")] = (1 .- uninsured_plans) .* 
			#	(households[household_indices,:perc_asian]);
			#X[:,findall(all_covariates .== "intercept_black")] = (1 .- uninsured_plans) .* 
			#	(households[household_indices,:perc_black]);	
			#X[:,findall(all_covariates .== "intercept_hispanic")] = (1 .- uninsured_plans) .* 
			#	(households[household_indices,:perc_hispanic]);
			#X[:,findall(all_covariates .== "intercept_other")] = (1 .- uninsured_plans) .* 
			#	(households[household_indices,:perc_other]);	
				
			# Premium interactions - Multiply by premium later
			X[:,findall(all_covariates .== "Premium_250to400")[1]] = (1 .- uninsured_plans) .* 
				((households[household_indices,:FPL] .> 2.5) .& (households[household_indices,:FPL] .<= 4)) ;	
			X[:,findall(all_covariates .== "Premium_gt400")[1]] = (1 .- uninsured_plans) .* 
				(households[household_indices,:FPL] .> 4) ; 
			X[:,findall(all_covariates .== "Premium_male")[1]] = (1 .- uninsured_plans) .*
				households[household_indices,:perc_male];
			X[:,findall(all_covariates .== "Premium_0to34")[1]] = (1 .- uninsured_plans) .* 
				(households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] + households[household_indices,:perc_26to34]);	
			X[:,findall(all_covariates .== "Premium_35to54")[1]] = (1 .- uninsured_plans) .* 
				(households[household_indices,:perc_35to44] + households[household_indices,:perc_45to54]);
			X[:,findall(all_covariates .== "Premium_family")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:household_size] .> 1);
			X[:,findall(all_covariates .== "Premium_Asian")[1]] = (1 .- uninsured_plans) .* 
				(households[household_indices,:perc_asian]);
			X[:,findall(all_covariates .== "Premium_Black")[1]] = (1 .- uninsured_plans) .* 
				(households[household_indices,:perc_black]);	
			X[:,findall(all_covariates .== "Premium_Hispanic")[1]] = (1 .- uninsured_plans) .* 
				(households[household_indices,:perc_hispanic]);
			X[:,findall(all_covariates .== "Premium_Other")[1]] = (1 .- uninsured_plans) .* 
				(households[household_indices,:perc_other]);	
			
			if add_complex_interactions
				X[:,findall(all_covariates .== "Premium_250to400_0to34")[1]] = (1 .- uninsured_plans) .* ((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) .*
					(households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]); 
				X[:,findall(all_covariates .== "Premium_gt400_0to34")[1]] = (1 .- uninsured_plans) .* (households[household_indices,:FPL] .> 4) .*
					(households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]);  
				X[:,findall(all_covariates .== "Premium_250to400_35to54")[1]] = (1 .- uninsured_plans) .* ((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) .*
					(households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]); 
				X[:,findall(all_covariates .== "Premium_gt400_35to54")[1]] = (1 .- uninsured_plans) .* (households[household_indices,:FPL] .> 4) .*
					(households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]);  
			end
			
			# AV interactions
			if add_av_interactions
				X[:,findall(all_covariates .== "AV_250to400")[1]] = X[:,av_ind] .* ((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) ; 
				X[:,findall(all_covariates .== "AV_gt400")[1]] = X[:,av_ind] .* (households[household_indices,:FPL] .> 4); 
				X[:,findall(all_covariates .== "AV_0to17")[1]] = X[:,av_ind] .* (households[household_indices,:perc_0to17]);
				X[:,findall(all_covariates .== "AV_18to34")[1]] = X[:,av_ind] .* (households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]);
				X[:,findall(all_covariates .== "AV_35to54")[1]] =  X[:,av_ind] .* (households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]);
				X[:,findall(all_covariates .== "AV_male")[1]] = X[:,av_ind] .* (households[household_indices,:perc_male]);
				X[:,findall(all_covariates .== "AV_family")[1]] = X[:,av_ind] .* (households[household_indices,:household_size] .> 1);
				X[:,findall(all_covariates .== "AV_Asian")[1]] = X[:,av_ind] .* (households[household_indices,:perc_asian]);
				X[:,findall(all_covariates .== "AV_Black")[1]] = X[:,av_ind] .* (households[household_indices,:perc_black]);
				X[:,findall(all_covariates .== "AV_Hispanic")[1]] = X[:,av_ind] .* (households[household_indices,:perc_hispanic]);
				X[:,findall(all_covariates .== "AV_Other")[1]] = X[:,av_ind] .* (households[household_indices,:perc_other]);
				
				if add_complex_interactions
					X[:,findall(all_covariates .== "AV_250to400_0to34")[1]] = X[:,av_ind] .* ((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) .*
						(households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]); 
					X[:,findall(all_covariates .== "AV_gt400_0to34")[1]] = X[:,av_ind] .* (households[household_indices,:FPL] .> 4) .*
						(households[household_indices,:perc_0to17] .+ households[household_indices,:perc_18to25] .+ households[household_indices,:perc_26to34]);  
					X[:,findall(all_covariates .== "AV_250to400_35to54")[1]] = X[:,av_ind] .* ((households[household_indices,:FPL] .> 2.50) .& (households[household_indices,:FPL] .<= 4)) .*
						(households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]); 
					X[:,findall(all_covariates .== "AV_gt400_35to54")[1]] = X[:,av_ind] .* (households[household_indices,:FPL] .> 4) .*
						(households[household_indices,:perc_35to44] .+ households[household_indices,:perc_45to54]);  
				end
			end
	
			
			# Market fixed effects
			X[:,findall(all_covariates .== "intercept_ra2")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 2);
			X[:,findall(all_covariates .== "intercept_ra3")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 3);
			X[:,findall(all_covariates .== "intercept_ra4")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 4);
			X[:,findall(all_covariates .== "intercept_ra5")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 5);
			X[:,findall(all_covariates .== "intercept_ra6")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 6);
			X[:,findall(all_covariates .== "intercept_ra7")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 7);
			X[:,findall(all_covariates .== "intercept_ra8")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 8);
			X[:,findall(all_covariates .== "intercept_ra9")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 9);
			X[:,findall(all_covariates .== "intercept_ra10")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 10);
			X[:,findall(all_covariates .== "intercept_ra11")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 11);
			X[:,findall(all_covariates .== "intercept_ra12")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 12);
			X[:,findall(all_covariates .== "intercept_ra14")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 14);
			X[:,findall(all_covariates .== "intercept_ra15")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 15);
			X[:,findall(all_covariates .== "intercept_ra16")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 16);
			X[:,findall(all_covariates .== "intercept_ra17")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 17);
			X[:,findall(all_covariates .== "intercept_ra18")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 18);
			X[:,findall(all_covariates .== "intercept_ra19")[1]] = (1 .- uninsured_plans) .*
				(households[household_indices,:rating_area] .== 19);
			
			# Insurer-market fixed effects
			
			if add_insurer_market_fe
			
				X[:,findall(all_covariates .== "intercept_ra2_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 2) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra2_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 2) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra2_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 2) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];	
				X[:,findall(all_covariates .== "intercept_ra2_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 2) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];

				X[:,findall(all_covariates .== "intercept_ra3_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 3) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra3_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 3) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra3_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 3) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];

				X[:,findall(all_covariates .== "intercept_ra4_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 4) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra4_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 4) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra4_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 4) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];	
				X[:,findall(all_covariates .== "intercept_ra4_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 4) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];

				X[:,findall(all_covariates .== "intercept_ra5_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 5) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra5_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 5) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra5_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 5) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];	
						
				X[:,findall(all_covariates .== "intercept_ra6_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 6) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra6_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 6) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
							
				X[:,findall(all_covariates .== "intercept_ra7_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 7) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra7_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 7) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra7_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 7) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];	
				X[:,findall(all_covariates .== "intercept_ra7_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 7) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];	
							
				X[:,findall(all_covariates .== "intercept_ra8_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 8) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra8_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 8) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra8_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 8) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];	
				X[:,findall(all_covariates .== "intercept_ra8_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 8) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];
						
				X[:,findall(all_covariates .== "intercept_ra9_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 9) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra9_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 9) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
						
				X[:,findall(all_covariates .== "intercept_ra10_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 10) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra10_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 10) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra10_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 10) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];	
						
				X[:,findall(all_covariates .== "intercept_ra11_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 11) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra11_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 11) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
						
				X[:,findall(all_covariates .== "intercept_ra12_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 12) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra12_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 12) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
						
				X[:,findall(all_covariates .== "intercept_ra14_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 14) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra14_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 14) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra14_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 14) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];
						
				X[:,findall(all_covariates .== "intercept_ra15_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 15) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra15_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 15) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra15_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 15) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];
				X[:,findall(all_covariates .== "intercept_ra15_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 15) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];
						
				X[:,findall(all_covariates .== "intercept_ra16_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 16) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra16_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 16) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];							
				X[:,findall(all_covariates .== "intercept_ra16_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 16) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];
				X[:,findall(all_covariates .== "intercept_ra16_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 16) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];
						
				X[:,findall(all_covariates .== "intercept_ra17_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 17) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra17_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 17) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra17_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 17) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];
				X[:,findall(all_covariates .== "intercept_ra17_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 17) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];
					
				X[:,findall(all_covariates .== "intercept_ra18_Anthem")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 18) .*
					X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra18_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 18) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra18_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 18) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];
				X[:,findall(all_covariates .== "intercept_ra18_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 18) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];	
							
				#X[:,findall(all_covariates .== "intercept_ra19_Anthem")[1]] =
				#	(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 19) .*
				#	X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "intercept_ra19_Blue_Shield")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 19) .*
					X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "intercept_ra19_Kaiser")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 19) .*
					X[:,findall(all_covariates .== "Kaiser")[1]];
				X[:,findall(all_covariates .== "intercept_ra19_Health_Net")[1]] =
					(1 .- uninsured_plans) .* (households[household_indices,:rating_area] .== 19) .*
					X[:,findall(all_covariates .== "Health_Net")[1]];
			end
			
			# Residual from reduced form stage
			if control_function
				X[:,findall(all_covariates .== "residual")[1]] = convert(Array,data[!,:residual]);
			end
			
			if random_parameters
				if !inattention
					X[:,findall(all_covariates .== "intercept_random")[1]] = X[:,findall(all_covariates .== "intercept")[1]];
				end
				X[:,findall(all_covariates .== "Anthem_random")[1]] = X[:,findall(all_covariates .== "Anthem")[1]];
				X[:,findall(all_covariates .== "Blue_Shield_random")[1]] = X[:,findall(all_covariates .== "Blue_Shield")[1]];
				X[:,findall(all_covariates .== "Kaiser_random")[1]] = X[:,findall(all_covariates .== "Kaiser")[1]];
				X[:,findall(all_covariates .== "Health_Net_random")[1]] = X[:,findall(all_covariates .== "Health_Net")[1]];
				if !add_av_interactions
					X[:,findall(all_covariates .== "AV_random")[1]] = X[:,findall(all_covariates .== "AV")[1]];
				end
			end
			
			
		return(X)
	end
	
	# Form X covariate matrix
	function form_Xsearch_matrix(search_covariates,households)	

		Xsearch = zeros(Float64,size(households)[1],length(search_covariates));
		Xsearch[:,findall(search_covariates .== "premium_shock")[1]] = households[!,:premium_shock];
		Xsearch[:,findall(search_covariates .== "hh_250to400")[1]] = ((households[!,:FPL] .> 2.50) .& (households[!,:FPL] .<= 4)) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_gt400")[1]] = (households[!,:FPL] .> 4) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_0to17")[1]] = (households[!,:perc_0to17]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_18to34")[1]] = (households[!,:perc_18to25] .+ households[!,:perc_26to34]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_35to54")[1]] = (households[!,:perc_35to44] .+ households[!,:perc_45to54]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Black")[1]] = (households[!,:perc_black]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Hispanic")[1]] = (households[!,:perc_hispanic]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Asian")[1]] = (households[!,:perc_asian]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Other")[1]] = (households[!,:perc_other]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_male")[1]] = (households[!,:perc_male]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_family")[1]] = (households[!,:household_size] .> 1) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2015")[1]] = (households[!,:year] .== 2015) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2016")[1]] = (households[!,:year] .== 2016) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2017")[1]] = (households[!,:year] .== 2017) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2018")[1]] = (households[!,:year] .== 2018) .* 1; 
		Xsearch[:,findall(search_covariates .== "intercept_search")[1]] .= 1; 
		
		return(Xsearch)
	end
	
	# Covariate Matrix for Claims Moments (data)
	
	function form_X_claims_matrix(moments_data,claims_moments_covariates)
		
		X_claims_data = zeros(size(moments_data)[1],length(claims_moments_covariates));	
		X_claims_data[:,findall(claims_moments_covariates .== "intercept")[1]] .= 1;
		X_claims_data[:,findall(claims_moments_covariates .== "HMO")[1]] = moments_data[!,:HMO];
		X_claims_data[:,findall(claims_moments_covariates .== "log_risk_score")[1]] = moments_data[!,:log_risk_score];
		X_claims_data[:,findall(claims_moments_covariates .== "trend")[1]] = moments_data[!,:trend];
		X_claims_data[:,findall(claims_moments_covariates .== "Anthem")[1]] = moments_data[!,:Anthem];
		X_claims_data[:,findall(claims_moments_covariates .== "Blue_Shield")[1]] = moments_data[!,:Blue_Shield];
		X_claims_data[:,findall(claims_moments_covariates .== "Kaiser")[1]] = moments_data[!,:Kaiser];
		X_claims_data[:,findall(claims_moments_covariates .== "Health_Net")[1]] = moments_data[!,:Health_Net];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra2")[1]] = moments_data[!,:share_ra2];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra3")[1]] = moments_data[!,:share_ra3];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra4")[1]] = moments_data[!,:share_ra4];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra5")[1]] = moments_data[!,:share_ra5];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra6")[1]] = moments_data[!,:share_ra6];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra7")[1]] = moments_data[!,:share_ra7];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra8")[1]] = moments_data[!,:share_ra8];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra9")[1]] = moments_data[!,:share_ra9];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra10")[1]] = moments_data[!,:share_ra10];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra11")[1]] = moments_data[!,:share_ra11];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra12")[1]] = moments_data[!,:share_ra12];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra13")[1]] = moments_data[!,:share_ra13];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra14")[1]] = moments_data[!,:share_ra14];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra15")[1]] = moments_data[!,:share_ra15];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra16")[1]] = moments_data[!,:share_ra16];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra17")[1]] = moments_data[!,:share_ra17];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra18")[1]] = moments_data[!,:share_ra18];
		X_claims_data[:,findall(claims_moments_covariates .== "share_ra19")[1]] = moments_data[!,:share_ra19];	
		
		return(X_claims_data)
	end
	
	# Covariate Matrix for Claims Moments (only exchange data)
	
	function form_X_CLMS_matrix(plans_pmt,hdata,claims_moments_covariates,theta_RS)
	
		X_CLMS = zeros(size(plans_pmt)[1],length(claims_moments_covariates));
		X_CLMS[:,findall(claims_moments_covariates .== "intercept")[1]] .= 1;
		X_CLMS[:,findall(claims_moments_covariates .== "log_risk_score")[1]] = X_RS * theta_RS;
		X_CLMS[:,findall(claims_moments_covariates .== "HMO")[1]] = plans_pmt[!,:HMO];
		X_CLMS[:,findall(claims_moments_covariates .== "trend")[1]] = plans_pmt[!,:year] .- 2014;
		X_CLMS[:,findall(claims_moments_covariates .== "Anthem")[1]] = plans_pmt[!,:Anthem];
		X_CLMS[:,findall(claims_moments_covariates .== "Blue_Shield")[1]] = plans_pmt[!,:Blue_Shield];
		X_CLMS[:,findall(claims_moments_covariates .== "Kaiser")[1]] = plans_pmt[!,:Kaiser];
		X_CLMS[:,findall(claims_moments_covariates .== "Health_Net")[1]] = plans_pmt[!,:Health_Net];
			
		# Add market fixed effect
		markets = plans_pmt[!,:rating_area];		
		for m in 1:length(all_mkts)
			X_CLMS[markets .== all_mkts[m],length(claims_moments_covariates) - length(rating_area_variables) + m] .= 1; 
		end	
		
		return(X_CLMS)
	end
	
	# Covariate Matrix for Risk Score Moments (only exchange data)
	
	function create_risk_score_share_variable_matrix(X,households)
	
		rs_share_variable_values = zeros(size(hdata)[1],length(risk_score_share_variables));
		
		if all_rs_variables
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]] = 
				(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_35to54")[1]] = 
				(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_250to400")[1]] = 
				((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_gt400")[1]] = 
				(households[household_numbers,:FPL] .> 4) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Asian")[1]] = 
				households[household_numbers,:perc_asian] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Black")[1]] = 
				households[household_numbers,:perc_black] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]] = 
				households[household_numbers,:perc_hispanic] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Other_Race")[1]] = 
				households[household_numbers,:perc_other] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]] = 
				households[household_numbers,:perc_male] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]] = 
				(households[household_numbers,:household_size] .> 1) .* 1;
		elseif int_rs_variables
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]] = 
				(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_minority")[1]] = 
				(households[household_numbers,:perc_asian] .+ households[household_numbers,:perc_black] .+ 
				households[household_numbers,:perc_hispanic] .+ households[household_numbers,:perc_other]) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]] = 
				households[household_numbers,:perc_male] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]] = 
				(households[household_numbers,:household_size] .> 1) .* 1;
		elseif hisp_flag
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]] = 
				(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]] = 
				households[household_numbers,:perc_hispanic] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]] = 
				households[household_numbers,:perc_male] .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_family")[1]] = 
				(households[household_numbers,:household_size] .> 1) .* 1;
		elseif agegender_flag
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_0to34")[1]] = 
				(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_male")[1]] = 
				households[household_numbers,:perc_male] .* 1;
		elseif old_flag
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_18to54")[1]] = 
				(households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34] .+ 
					households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]) .* 1;
			rs_share_variable_values[:,findall(risk_score_share_variables .== "share_Hispanic")[1]] = 
				households[household_numbers,:perc_hispanic] .* 1;
		end
		
		return(rs_share_variable_values)
	end
	
	function form_X_RS_matrix(plans_pmt,hdata,household_numbers,risk_score_covariates)
	
		X_RS = zeros(size(plans_pmt)[1],length(risk_score_covariates));
		X_RS[:,findall(risk_score_covariates .== "intercept")[1]] .= 1;
		if av_rs_flag
			X_RS[:,findall(risk_score_covariates .== "AV")[1]] = plans_pmt[!,:AV];
		elseif !separate_metal_flag
			X_RS[:,findall(risk_score_covariates .== "Silver")[1]] = (plans_pmt[!,:AV] .== 0.7) .* 1;
			X_RS[:,findall(risk_score_covariates .== "Gold")[1]] = (plans_pmt[!,:AV] .== 0.8) .* 1;
			X_RS[:,findall(risk_score_covariates .== "Platinum")[1]] = (plans_pmt[!,:AV] .== 0.9) .* 1;
		end
		
		if inattention
			probability_weights = hdata[:,findall(hdata_fields .== "real_prob")[1]] .* hdata[:,findall(hdata_fields .== "weight")[1]];
		else
			probability_weights = hdata[:,findall(hdata_fields .== "real_instrumented_prob")[1]] .* hdata[:,findall(hdata_fields .== "weight")[1]];
		end
	
		if all_rs_variables
			weight_object = DataFrame(plan_name = data_pmt_names,
							probability_weights = probability_weights,
							probability_weights_0to34 = probability_weights .* (households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]),
							probability_weights_35to54 = probability_weights .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]),
							probability_weights_250to400 = probability_weights .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)),
							probability_weights_gt400 = probability_weights .* (households[household_numbers,:FPL] .> 4),
							probability_weights_Male = probability_weights .* households[household_numbers,:perc_male],
							probability_weights_Family = probability_weights .* (households[household_numbers,:household_size] .> 1),
							probability_weights_Asian = probability_weights .* households[household_numbers,:perc_asian],
							probability_weights_Black = probability_weights .* households[household_numbers,:perc_black],
							probability_weights_Hispanic = probability_weights .* households[household_numbers,:perc_hispanic],
							probability_weights_Other = probability_weights .* households[household_numbers,:perc_other]);
		elseif int_rs_variables
			weight_object = DataFrame(plan_name = data_pmt_names,
							probability_weights = probability_weights,
							probability_weights_0to34 = probability_weights .* (households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]),
							probability_weights_Male = probability_weights .* households[household_numbers,:perc_male],
							probability_weights_Family = probability_weights .* (households[household_numbers,:household_size] .> 1),
							probability_weights_Minority = probability_weights .* 
								(households[household_numbers,:perc_asian] .+ households[household_numbers,:perc_black] .+ 
									households[household_numbers,:perc_hispanic] .+ households[household_numbers,:perc_other]));
		elseif hisp_flag
			weight_object = DataFrame(plan_name = data_pmt_names,
							probability_weights = probability_weights,
							probability_weights_0to34 = probability_weights .* (households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]),
							probability_weights_Male = probability_weights .* households[household_numbers,:perc_male],
							probability_weights_Family = probability_weights .* (households[household_numbers,:household_size] .> 1),
							probability_weights_Hispanic = probability_weights .* households[household_numbers,:perc_hispanic]);
		elseif agegender_flag
			weight_object = DataFrame(plan_name = data_pmt_names,
							probability_weights = probability_weights,
							probability_weights_0to34 = probability_weights .* (households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]),
							probability_weights_Male = probability_weights .* households[household_numbers,:perc_male]);
		elseif old_flag
			weight_object = DataFrame(plan_name = data_pmt_names,
							probability_weights = probability_weights,
							probability_weights_18to54 = probability_weights .* (households[household_numbers,:perc_18to25] .+ 
								households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]),
							probability_weights_Hispanic = probability_weights .* households[household_numbers,:perc_hispanic]);
		end
		
		total_plan_weights = sort!(combine(groupby(weight_object,:plan_name), :probability_weights => sum),[:plan_name])[!,:probability_weights_sum];
		
		if all_rs_variables
			X_RS[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_0to34 => sum),[:plan_name])[!,:probability_weights_0to34_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_35to54")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_35to54 => sum),[:plan_name])[!,:probability_weights_35to54_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_250to400")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_250to400 => sum),[:plan_name])[!,:probability_weights_250to400_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_gt400")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_gt400 => sum),[:plan_name])[!,:probability_weights_gt400_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Male => sum),[:plan_name])[!,:probability_weights_Male_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Family => sum),[:plan_name])[!,:probability_weights_Family_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_Asian")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Asian => sum),[:plan_name])[!,:probability_weights_Asian_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_Black")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Black => sum),[:plan_name])[!,:probability_weights_Black_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Hispanic => sum),[:plan_name])[!,:probability_weights_Hispanic_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_Other_Race")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Other => sum),[:plan_name])[!,:probability_weights_Other_sum] ./ total_plan_weights;
		elseif int_rs_variables
			X_RS[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_0to34 => sum),[:plan_name])[!,:probability_weights_0to34_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Male => sum),[:plan_name])[!,:probability_weights_Male_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Family => sum),[:plan_name])[!,:probability_weights_Family_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_minority")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Minority => sum),[:plan_name])[!,:probability_weights_Minority_sum] ./ total_plan_weights;
		elseif hisp_flag
			X_RS[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_0to34 => sum),[:plan_name])[!,:probability_weights_0to34_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Male => sum),[:plan_name])[!,:probability_weights_Male_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_family")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Family => sum),[:plan_name])[!,:probability_weights_Family_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Hispanic => sum),[:plan_name])[!,:probability_weights_Hispanic_sum] ./ total_plan_weights;
		elseif agegender_flag
			X_RS[:,findall(risk_score_covariates .== "share_0to34")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_0to34 => sum),[:plan_name])[!,:probability_weights_0to34_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_male")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Male => sum),[:plan_name])[!,:probability_weights_Male_sum] ./ total_plan_weights;
		elseif old_flag
			X_RS[:,findall(risk_score_covariates .== "share_18to54")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_18to54 => sum),[:plan_name])[!,:probability_weights_18to54_sum] ./ total_plan_weights;
			X_RS[:,findall(risk_score_covariates .== "share_Hispanic")[1]] = 
				sort!(combine(groupby(weight_object,:plan_name), :probability_weights_Hispanic => sum),[:plan_name])[!,:probability_weights_Hispanic_sum] ./ total_plan_weights;
		end
		
		return(X_RS)
	end
	
	